%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                        Employment Distribution                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Dy1 = Employ_Dist(p, F1, AgeDensity_N, Deltagrid, LPNwMLam)

% --------------------------- Description --------------------------------- 

    % This function calculates the distribution of firm size (in terms of
    % employees) given the firm's number of products.

    % For an overview of this funciton, please see Section (RP-2.5)

%--------------------------------------------------------------------------
%              1. Conditional Distribution with Competitors               %
%--------------------------------------------------------------------------

    % Define the number of points in the gap distribution
    Deltagrid_pts = length(Deltagrid);

    % We denote the distribution with competitors by Fc
    Fc = F1;

    % We take differences to obtain the associated density
    Fc_dens = zeros(Deltagrid_pts, p.n3);
    Fc_dens(:, 1:end-1) = diff(Fc, 1, 2);
    Fc_dens = [diff(Fc_dens, 1, 1); zeros(1, p.n3)];

    % Compute the marginal distribution - denominator in Equation (RP-39)
    Fdelta_marginal = sum(Fc_dens, 2);

    % Compute the conditional distribution accordingly - see Equation (RP-39)
    Fc_cond = cumsum(Fc_dens, 2) ./ repmat(Fdelta_marginal, 1, p.n3);

    % For the last possible gap, there are NaNs - repeat previous line
    Fc_cond(end, :) = Fc_cond(end-1, :);
    
%--------------------------------------------------------------------------
%            2. Conditional Distribution without Competitors              %
%--------------------------------------------------------------------------

    % Compute the implied step size in the productivity grid
    delta_q = p.qgrid(2) - p.qgrid(1);

    % Create a grid for firm age that matches the one used for the conditional distribution
    TimePeriods = p.maxTime/p.periodlength + 1;
    avec = linspace(0, p.maxTime, TimePeriods);

    % Compute the step implied in the firm age grid
    delta_a = avec(2) - avec(1);

    % Preallocate a matrix for the conditional distribution F^{NC}(q) over possible firm age
    FNC_a = zeros(TimePeriods, p.n3);

    % Transformed productivity equivalent to minimum relative efficiency
    MinEff = log(p.wmin);

    % Truncate the grid for efficiency to set cdf to zero before MinEff
    grid_aux = max(MinEff, p.qgrid);

    % Initial condition for the differential equation - see Equation (RP-36)
    FNC_a(1, :) = 1 - (p.wmin*exp(-grid_aux)).^p.varrho;

    % Over each value of firm age in the grid...
    for a = 2:TimePeriods

        % ... compute the derivative wrt q over the previous line - see Equation (RP-37)
        dervec = diff(FNC_a(a-1,:));
        dfdq = [dervec(1), dervec]/delta_q;

        % ... compute the derivative with respect to firm age - see Equation (RP-37)
        dfda = -dfdq*(p.In-p.gQ) - FNC_a(a-1,:)*(p.tau+p.d0);

        % ... computes the next row of the distribution  - see Equation (RP-37)
        FNC_a(a,:) = FNC_a(a-1,:) + dfda*delta_a;

    end

    % Rescale to get the conditional distribution - see Equation (RP-38)
    FNCond_a = FNC_a ./ FNC_a(:, end); 

%--------------------------------------------------------------------------
%                     3. Defining the appropriate grids                   %
%--------------------------------------------------------------------------

    % Varieties without competitors: grid for log(y) implied in our q grid
    ygrid_endog_NC = (p.s-1)*p.qgrid_adj - p.s*log(p.s/(p.s-1)) + log(LPNwMLam);
    ygrid_endog_NC = repmat(ygrid_endog_NC, TimePeriods, 1);

    % Varieties with competitors: grid for log(y) implied in our q grid
    GapGrid = p.l*exp(p.In*avec);
    MarkupGrid = min(p.s/(p.s-1), GapGrid);
    ygrid_endog_C = repmat((p.s-1)*p.qgrid_adj,TimePeriods,1) - p.s*log(repmat(MarkupGrid', 1, p.n3))+ log(LPNwMLam);

    % The target (and exogenous) grid for employment comes from Params
    yvec = p.yvec;
    y_grid = repmat(yvec, TimePeriods, 1);

    % Create matrices to store the interpolated distributions
    FNC_cond_grid = zeros(TimePeriods, p.ypoints);
    FC_cond_grid = zeros(TimePeriods, p.ypoints);

    % Determine the mesh grids where Fc_cond is defined
    [Mesh1_start, Mesh2_start] = meshgrid(Deltagrid, p.qgrid_adj);

    % Determine the mesh grids where we want Fc_cond
    [Mesh1_end, Mesh2_end] = meshgrid(GapGrid, p.qgrid_adj);

    % Interpolate Fc_cond to the desired grid
    Fc_interp = interp2(Mesh1_start, Mesh2_start, Fc_cond', Mesh1_end, Mesh2_end, 'spline', 0);
    Fc_interp = Fc_interp';
    clear Mesh1_start Mesh2_end Mesh1_start Mesh2_end 

    % Interpolate the conditional distribution we derived to the grid
    for A = 1:TimePeriods

        % Interpolate the distribution without competitors to the grid
        FNC_cond_grid(A,:) = (1-p.alpha)*interp1(ygrid_endog_NC(A,:), FNCond_a(A,:), y_grid(A,:));

        % Interpolate the distribution with competitors to the grid
        FC_cond_grid(A,:) = p.alpha*interp1(ygrid_endog_C(A,:), Fc_interp(A,:), y_grid(A,:));    

    end

    % Since ygrid starts before the endogenous grids, there are NaNs
    FNC_cond_grid(isnan(FNC_cond_grid)) = 0;
    FC_cond_grid(isnan(FC_cond_grid)) = 0;

%--------------------------------------------------------------------------
%                  4. Computing Product-Level Employment                  %
%--------------------------------------------------------------------------

    % Preallocate a matrix to store the D1_{n}(y) function
    Dy1 = zeros(p.ypoints, p.maxprods);

    % For every number of products a firm might have...
    for n = 1:p.maxprods

        % ... we redimension the product age distribution for calculation - see Equations (RP-34) and (RP-35)
        Density_aux = repmat(AgeDensity_N(:,n), 1, p.ypoints);

        % ... we compute the integrand for all possible product age a - see Equations (RP-34) and (RP-35)
        int_term = Density_aux .* (FC_cond_grid + FNC_cond_grid);

        % ... sum over product age to obtain D_n^1(y) - see Equations (RP-34) and (RP-35)
        Dy1(:, n) = sum(int_term, 1);

    end

    % Since there is interpolation out the imp_grids, there are NaNs
    Dy1(isnan(Dy1)) = 0;

end


   

   